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On large-scales, comparable to the horizon, the observable clustering properties of galaxies are affected by 
various general relativistic effects. To calculate these effects one needs to consistently solve for the metric, 
densities and velocities in a specific coordinate system or gauge. The method of choice for simulating large- 
scale structure is numerical N-body simulations which are performed in the Newtonian limit. Even though one 
might worry that the use of the Newtonian approximation would make it impossible to use these simulations to 
compute properties on very large-scales, we show that the simulations are still solving the dynamics correctly 
even for long modes and we give formulas to obtain the position of particles in the conformal Newtonian gauge 
given the positions computed in the simulation. We also give formulas to convert from the output coordinates 
of N-body simulations to the observable coordinates of the particles. 
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■ The study of the fluctuations in the distribution of matter in the 
^> Universe and its evolution through cosmic history has become 
r \one of the major tools in cosmology. The properties and time 
dvolution of the large-scale structure depend on the cosmologi- 
cal parameters and on the initial conditions for the hot big bang. 
OMany of the parameters of the currently favored cosmological 
Q model have been determined by matching the observed proper- 
?— I ties of the distribution of mass through cosmic history with the 
c/) model calculations. 

^ ^ ^ Galaxies serve as tracers of the underlying matter distribution. 

Significant efforts have been made to understand their connec- 
(T^tion iQl-Q] and to generate estimates for cosmological parame- 
J>ters from the recovered matter power spectrum |5]. In the last 
decades, redshift surveys such as the Sloan Digital Sky Survey 
•^(SDSS) (@] and flie Two-degree Field Galaxy Redshift Survey 
XJ(2dFGRS) fh have resulted in detailed maps of the large-scale 
• distribution of galaxies across very large volumes. The future 
promises even larger surveys as a result of efforts to improve 
measures of the so-called baryon acoustic oscillation signal in 
I the clustering of matter to further constrain the properties of the 
• .dark energy |l8l[9[]. Surveys are beginning to probe scales com- 
parable to the horizon at the redshift of the galaxies being ob- 
k> served. 

^ Until recently, predictions for observables in galaxy surveys 
5^ had been done entirely in the Newtonian limit. Matsubara iflOll 
included gravitational lensing effects on the correlation func- 
tions of galaxies and quasars as applied to SDSS. More recently, 
work by Yoo et al. 11 1.1 . Yoo LlZl made a more detailed treat- 



ment of general relativistic effects. These become relevant as 
the scales probed by the survey approach the horizon scale. The 
overdensity in the galaxy distribution Jobs is given by 
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where H is the Hubble constant; 5z, Sr and SVl are the fluc- 
tuations in the redshift, distance along the line of sight, and lu- 
minosity distance relative to the unperturbed universe; k is the 
lensing convergence; p gives the slope of the galaxy luminos- 
ity function; b is the bias; the direction of propagation of the 
photon; and A, Bi, D and Eij are metric components: 



ds' = 



a' [(14 



2 A) dif -2a^Bi dr] dx' 



(2) 



2D)g,^ 



2Eii\ dx'dx' 



with gij the metric tensor for three-space in a homogeneous 
universe and ?/ = J dt/a{t) is the conformal time in terms of 
the scale factor. These formulas exhibit many of the relativis- 
tic effects that are common in calculation of the anisotropics in 
the cosmic microwave background (CMB). For example the ob- 
served redshift of a source is given by 
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where the prime indicates the derivative with respect to confor- 
mal time, the vertical bar is the covariant derivative with respect 
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to gij, Ts — r{zs) is the comoving line of sight distance to the 
source galaxies at Zg, u^e* is the line of sight peculiar veloc- 
ity, and flo and are the values of the scale factor at the time 
of observation and light-emission respectively. The first square 
bracket represents the redshift-space distortion by peculiar ve- 
locities, frame dragging, and gravitational redshift, respectively. 
The first round bracket in the integral also represents the gravi- 
tational redshift, arising from the net difference in gravitational 
potential due to its time evolution for the duration of photon 
propagation, and this effect is referred to as the integrated Sachs- 
Wolfe (ISW) effect in the CMB literature. The last terms in the 
integral represent the tidal effect from the frame dragging and 
the ISW effect from the time evolution of the primordial gravity 
waves. 

The complete set of formulas needed to predict the observed 
clustering properties of galaxies on very large-scales can be 
found in Yoo et al. [11]. It is clear that the calculation requires 
consistently solving the general relativity dynamics in a partic- 
ular coordinate system or gauge. On the other hand N-body 
simulations are the method of choice to compute predictions for 
the large-scale distribution of galaxies but these simulations are 
done in the Newtonian limit. It is appropriate to ask how the 
output of simulations can be used to compute the different terms 
in Eq. ([U and even whether this can be done at all given that the 
simulations are run using Newtonian dynamics. 

The drive on the observational side to map larger and larger 
volumes of the Universe and the exponential increase in com- 
puter power have also resulted in numerical simulations of 
ever increasing size. Typical cosmological simulations evolve 
the particles starting at z ~ 100, when the size of the hori- 
zon is ^ l.SGpc. Box sizes vary and can be as large as 
^ 0.5 — 3Gpc comoving and the number of particles is of order 
^ 10^ — 10^*^. Examples of some of the biggest simulations to 
date are the Millennium Simulation 1 13] run in a box of comov- 
ing size 500/i^^Mpc, the Marenostrum Numerical Cosmology 
Project 1 14] and the Hubble Volume project 11151] run in a box 
of 3000 X 3000 X 30 Mpc. Further examples are found in Col- 
berg et al. iH], Park et al. |[l3l and recently in Cai et al. ifisll . 
Some of these simulations are started at an initial time when the 
horizon actually lies within the box. Clearly, we need a way to 
match cosmological N-body simulations with the general rela- 
tivistic variables in Eq. ([TJ. 

In addition to asking how to use the outputs of numerical sim- 
ulations to compute the various terms in Eq. ([T]i one may won- 
der if numerical simulations are solving the correct dynamical 
equations. We might suspect that the Newtonian simulations are 
working in the so-called conformal Newtonian gauge, in which 
the line element is given by 

ds^ = a^(?7)[-(l + '2.(t)N)drf + (1 - 2(j)N)5^jdx'dx^] (4) 

in the absence of anisotropic stress and where 4)^ coincides with 



the Newtonian potential only on small scales. In fact, the ana- 
logue of the Poisson equation in the conformal Newtonian gauge 
reads: 

V'^(j)N -in{(l)'N+n(l)N)=^T^Ga^5pN (5) 

and thus differs from the standard Poisson equation on large- 
scales. Here, T-L — a' /ais the conformal Hubble parameter. 

From Eq. (|5]l, it might appear that simulations are not solving 
correctly for the gravitational potential for scales comparable to 
or larger than the horizon. Previous work on this subject has fo- 
cused on comparing the general relativity equations to the New- 
tonian equations up to some given order in perturbation theory 
ifigll . We will show in this work that a more direct approach is 
possible. We will analyze the situation in detail and conclude 
that N-body simulations are solving for the potential correctly 
but that the location of the particles needs to be corrected if they 
are to be interpreted as the particle coordinates in the conformal 
Newtonian gauge. Finally we will give formulas to recover ob- 
servable coordinates directly in terms of the output of N-body 
simulations. 



II. EVOLUTION EQUATIONS 

As structure in the Universe develops the density contrast be- 
comes larger and larger exceeding unity at the so-called nonlin- 
ear scale. Properly modeling this process on small scales, of 
order the nonlinear scale or smaller requires numerical simula- 
tions. However, because the primordial curvature fluctuations, 
the seeds for structure formation, are so small, the nonlinear 
scale is significantly smaller than the horizon. As a result per- 
turbations in the space time remain very small, of order 10^^ or 
smaller'. 

Thus, to study structure formation we need only consider 
small perturbations to the Friedmann-Robertson- Walker metric 
and we can stay at linear order on those perturbations. In this pa- 
per we choose to work in the conformal Newtonian gauge with 
the line element given by: 

ds^ = a2(?/)[-(l + 2ijN)d7f + (1 - 2(l)N)5^jdx'dx^] (6) 

where ipM represents the Newtonian potential and (f)^, the New- 
tonian curvature. 



' On sufficiently small scales baryons can collapse to fomi relativistic objects 
such as neutron stars or black holes around which the space time metric fluc- 
tuations are large. This has negligible effects on the length scales considered 
in this paper 
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We stress that we are assuming that metric perturbations are 
small but we are not treating the density perturbations using per- 
turbation theory. The structure formation process also results in 
peculiar velocities for the particles. Because the nonlinear scale 
is well inside the horizon, these peculiar velocities are small, 
much smaller than the speed of light. In fact at the nonlinear 
scale peculiar velocities are of order the Hubble velocity for 
points separated by a distance of order the nonlinear scale. As 
a result the kinetic energies of particles do not source gravity in 
an appreciable way. 

Note that we are considering perturbations around the FRW 
metric so at lowest order the source for the gravitational poten- 
tials in 5p as opposed to the full p. The kinetic energy correc- 
tions are of order pv^. It is still true that pv^ <^ 6p on every 
scale. Thus it is safe to ignore the peculiar motions as a source 
of gravity. These terms are of course also neglected in numerical 
simulations run using the Newtonian approximation, but this is 
a negligible source of error. Including these terms is necessary 
if one wants to study the back-reaction of cosmological pertur- 
bations on to the expansion of the Universe [20], but we are not 
interested in this problem here. 

In the standard Newtonian approximation terms of order p(j) 
are also dropped as sources of gravity. This requires a little bit 
more thought in our case. Again at lowest order the source of 
gravity is 6p but this is no longer much larger than p(p on suffi- 
ciently large-scales. Thus we need to keep this term. However 
it is only p(f) that needs to be kept as of course Sp 3> Spcf) on all 
scales. We will now summarize the evolution equations under 
these approximations. 



A. Einstein equations 

In the conformal Newtonian gauge the Einstein equations 

Gf_i,y — SnGT^f, are reduced to: 



where ^ indicates derivatives with respect to the i coordinate 
and the background Friedmann equation for a flat Universe with 
cosmological constant A has been subtracted. 



— = -A.GT,+-, 



(11) 



B. Application to nonrelativistic particles 

As we mentioned when studying structure formation we are 
primarily interested in nonrelativistic matter Equation ( fTOb for 
i ^ i implies that the anisotropic stress is of order pv'^ and thus 
negligible in our approximation. As a result (pN = ipN- The 
gravitational potential satisfies: 
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The energy-momentum tensor for a set of cold dark matter par- 
ticles with mass ma [21] is given by: 
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where is the comoving velocity of the particles, dx'^'' /drj and 
g is the determinant of the metric. At linear order in the metric 
and to second order in the three-velocities, u*, this is given by 
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nent of the stress-energy tensor is related to the density of parti- 
cles and because metric perturbations are small, we can expand 
in powers of and remain at linear order. To order v^, the 00 
component is 

Tl^^-a-^'Y^mail + S^N+vDdDix^-Xam- (14) 
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As discussed we will neglect the term but need to keep the 
4)N as (fiMP is not negligible on large-scales. We obtain 

T° = -a-^(l + 30jv) "^aSoix: - Xa{v)) (15) 



AtiGo^ 



1 



hj^'^{4>N - i'N)] 



(9) 



(10) 



where only the 0jvP piece of the term proportional to (j)N ever 
makes any difference. Replacing in Eq. (fT2l i. 
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where p{x, rf) — a ^ niaSDix — Xa{ri)) is the density ob- 
tained by naively counting particles in cells at each time step. 
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Given the positions of the particles, Eq. ( fTSI l can be solved to 
obtain the Newtonian potential. 

The positions and velocities of the particles are advanced us- 
ing the geodesic equation. 



drj 



(17) 



Notice that the term 3(j)'j^dxa/dri is always negligible. 

m. INITIAL CONDITIONS 

In addition to the evolution equations we need to find the ini- 
tial conditions. This can be done at early enough time using 
linear theory. We define a growth function for the potential such 
that in the linear regime 



(t>N = b^{vi)(j)' 



N- 



(18) 



Given that the spatial components of the stress tensor for the 
dark matter are negligibly small, the gravitational potential in 
the matter era satisfies 



N 



0. 



(19) 



The solutions to this equation are well-known f^^], giving a 
constant and a decaying solution for the potential in the linear 
regime, which in terms of Eq. ( fTSl ) is 



9i 



(20) 



where Ci and C2 are constants. The decaying mode is abso- 
lutely negligible at the times of interest and without loss of gen- 
erality we choose Ci — 1. 

We can replace the right-hand side of the geodesic equation 
by the solution of the potential in the linear regime. 
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In a matter-dominated regime the homogeneous solutions to Eq. 
(ETJ are given by a constant vector and a decaying solution. 



Xh Bir]-^ + B2. 



(22) 
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the particular solution, we choose an ansatz 
— bs{r]) ijMxin), as it is usually done in the Zel'dovich 
approximation |23]. The labeling of 65(77) as such will become 
clear by the end of this section. The equation to be solved is. 
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The right-hand side of the last equation is independent of time, 
which implies that h'^ + T-ib'g — constant. The actual value 
of the constant is arbitrary, since it can be absorbed in ipi{xin)- 
Thus, bs oc f]^. (Adding a constant to this solution would not 
modify the subsequent steps of our paper and is only linked to 
the choice of initial time.) If we equate the factors that depend 
on the coordinates, ^pl{xin) oc —V(j)^j!}. 

To give the complete solution for the position of the particles 
we separate the constant vector B2 in two components, B2 = 
Xin + 6xin{xin), whcrc Xin are the positions of the particles if 
they started out distributed uniformly in a mesh. We discard the 
decaying term and give the position of the particles as a function 
of time in the linear regime. 
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To find the value of Sxin{xin), we resort to the Poisson equation 
(fTSI l. We can evolve p{x, rf) by means of a transformation of 
coordinates from the initial particle density 
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where 



dx 



is the Jacobian of the transformation and p is the 



initial uniform background density. The transformation of co- 
ordinates is given by Eq. ( |24] |. where 6x, bs{r]) and -01 are un- 
knowns. Since the perturbations are initially small, the density 
evolves as 
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We can define contributions to the density perturbation, 5, as 
related to the displacement fields by 
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The physical meaning of bs{ri) now becomes clear, as it is iden- 
tified with the growth function of density perturbations. Indeed, 
as we expect in the matter-dominated regime, we obtained that 
65(77) oc 0(77) oc ?7^. 

To determine we replace Eqns. (l26b . (l27b and (|28]) in Eq. 
(fTSI l to obtain: 



S^). (29) 



At the initial time the first term in the left-hand side cancels with 
the rightmost term in the right-hand side of the previous equa- 
tion. We can solve for (5"' from the remaining terms. 



N- 



(30) 
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In terms of their Fourier components, while (5™ oc (/i™ j,, the 
Zel'dovich component dependency is oc fc^H^^i/)^^ Well 
inside the horizon, in the limit kr] ^ 1, (5™ becomes negligible 
as compared to , the Newtonian density perturbation. When 
kr] 1 then ^ and as anticipated we cannot neglect 
this term. 

Finally, we can obtain 6xin by inverting Eq. (IZTl i in Fourier 
space. 



(Pk ik 
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In summary, to perform a cosmological simulation we need to 
evolve the position and velocities of particles using Eq. (fTTI ) and 
compute the potential using Eqs. Q and ( fTSl l. At the initial time 
the cold dark matter particles need to be displaced by an amount 
given in Eq. dSTT i. 



IV. COMPARISON TO NEWTONIAN COSMOLOGICAL 
SIMULATIONS 

Now that we have a consistent set of equations to solve we 
can compare them to those used in cosmological simulations to 
determine whether these simulations can be used to study very 
long wavelength modes or if they require some change. 

Cosmological Newtonian simulations solve for the potential 
by means of the Poisson equation 



(32) 



and move particles according to Newton's law expressed in co- 
moving coordinates. 
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This evolution equation is identical to the geodesic equation ( fTTI ) 
if the gravitational potential were computed correctly (as men- 
tioned before the term oc S^at' is negligible for modes both large 
and small compared to Hubble). Thus if the gravitational poten- 
tial is correct the particle positions are updated properly. 

It is important to determine if there are corrections to the grav- 
itational potential that become important on large-scales. The 
density that sources the Poisson equation in simulations is di- 
rectly computed by counting particles in cells, 

Psim ix,r/) ^ Msim ^ a'^^maSoix - Xail])). (34) 

a 

Even if the particle positions had been computed correctly, this 
"simulation density" differs from the density in the conformal 
Newtonian gauge by a factor (1 + S^at). 



Finally in standard cosmological simulation the particles are 
initially displaced making use of the Zel'dovich approximation, 
which in the Newtonian case takes the form. 



Xa = Xin + bs{r])'ljji{xin). 



(35) 



This differs from the displacements we calculated in the previous 
section - it is missing the Sxin- 

Thus at first sight it appears that the gravitational potential is 
not computed using the correct equation, that the density con- 
trast is missing a term and that the initial displacement of the 
particles is incorrect. We will now show that in fact all these dif- 
ferent "missing terms" cancel each other so that the gravitational 
potential is computed correctly. As a result, particle positions are 
also updated correctly. 

Let us look at the situation more carefully. For completeness 
let us also include a cosmological constant and start by restrict- 
ing ourselves to the linear regime as in any event the effects we 
are considering are only relevant on very large-scales. The rela- 
tivistic Poisson equation reads 



(36) 



where we have introduced the equation of state parameter co — 
p/p and we have written —V • Sxin = S"^ = SCin- Notice that 
3/2H^(l + u) = ATrGpdm with pdm the mean density for the 
matter. The terms in brackets on the right-hand side correspond 
to the density contrast calculated in the Newtonian simulations 
{5^) and the two missing corrections, the one proportional to 
and the one coming from the missing initial displacements (Cin)- 
It is useful to consider the comoving curvature defined as 



c 
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(37) 



It is well known that this comoving curvature remains constant in 
time on large-scales. For completeness we spell out the deriva- 
tion in the Appendix. Note that C,in defined above is nothing 
other than the initial value for this variable In the case of 
a perfect fluid C is constant on all scales larger than the sound 
horizon, or fc^c^ << TH? . In general what plays the role of 
is just the typical velocity dispersion that relates the magnitude 
of the spatial components of the energy-momentum tensor to the 
density. In our case, it is the velocity dispersion of the dark mat- 
ter particles induced by the growth of perturbations and thus 
remains constant all the way to the non linear scale. 

The relativistic corrections to the Newtonian Poisson equa- 
tion, V^0Ar — 3/2H^(l + oj)5^ , are evidenced from subtract- 
ing this expression from Eq. (l36T l. Notice that the difference 
between the relativistic terms in the right- and left-hand sides of 



6 



the Poisson equation is nothing other than: 

3H(0^ + n4>N) + ^n^i + u){4>N - Cm) = (38) 

Thus the additional terms in the relativistic Poisson equation 
cancel each other for modes larger than the nonlinear scale. Fur- 
thermore, because the nonlinear scale is well inside the horizon 
once the cancellation begins to fail the additional terms are very 
small, of order the ratio of the nonlinear scale to the horizon 
squared. Thus the gravitational potential in Newtonian simula- 
tions coincides with the one in the conformal Newtonian gauge 
even on very large-scales and thus the dark matter particles are 
being displaced correctly. 

We have showed that Newtonian N-body simulations calcu- 
late the correct gravitational potential and displace the particles 
correctly. The coordinates of the particles however are not the 
coordinates in the conformal Newtonian gauge as they are miss- 
ing the initial displacement. Thus the dictionary between con- 
formal Newtonian gauge variables and numerical simulations is: 

</>Af = 0sim, (39) 
m = Vsim, (40) 
XN = Xsim + Sxin- (41) 

Notice that the reason why the scheme worked was that the 
"sound horizon" scale, the scale out to which C was constant in 
time, was well inside the horizon. This scale is nothing other 
than the scale that dark matter particles can move since the time 
of the Big Bang as a result of the peculiar velocities they have. 
Because the dark matter particles are non-relativistic this dis- 
tance is well inside the horizon and the mistakes are negligible. 
Of course a simulation based on Newtonian physics could not 
work if particles are moving at an appreciable fraction of the 
speed of light. In general if there is an additional relativistic 
component the Newtonian simulations would not be computing 
things properly. The cosmological constant did not cause any 
problem because even though it is in some sense relativistic it 
is homogeneous so it does not contribute to the perturbation of 
the stress tensor. Thus as long as one is modeling nonrelativis- 
tic components or a relativistic component that does not cluster 
one is safe using the Newtonian approximation. Such a New- 
tonian approximation would not work, for example, during the 
radiation era where the effective sound speed of the dominant 
component of the energy density is very close to the speed of 
light. 



V. THE COMOVING GAUGE 

For completeness we will now show that Eq. ( fT2] l can be writ- 
ten making apparent gauge transformations between the comov- 
ing and the conformal Newtonian gauge. 

In Eq. ( l36b . the factor -^S^ is the density perturbation in 
Newtonian simulations. The expression in brackets is the con- 
formal Newtonian gauge density perturbation. Indeed rearrang- 
ing the terms we have 

\7^^N-5ni(b'N + nM = (42) 

As long as ( is constant, we can recognize in the right-hand side 
the gauge transformation for the density perturbation between 
the comoving and the conformal Newtonian gauge, 

SpN^5pc + S{^N-C)pi^+^)- (43) 

This is a well-known relation ll24ll and in this context it implies 
that the density perturbation that Newtonian simulations are ob- 
taining is the one in the comoving gauge. 

VI. OBSERVABLE COORDINATES 

In an inhomogeneous universe, the observed positions of the 
particles in the simulation are modified due to effects such as 
the Sachs-Wolfe effect, gravitational lensing, and peculiar ve- 
locities. The net result is that photons from a source follow a 
path that is perturbed with respect to the light cone of an ob- 
server in a homogeneous universe. Consider a comoving ob- 
server in an inhomogeneous universe with a velocity given by 
— {{1 — iI>m) / cLtV^ / a). The direction of observation is n, 
defined by {6, (p) in spherical coordinates, but due to the per- 
turbations to the photon path, the direction toward the point 
reached by the photon geodesic is actually s, corresponding to 
{9 + S9,(p + 6(p). We now take a particle with coordinates Xa{ri) 
(already taking into account the correction Sxin) and we want to 
know where it crosses the path of the photons going toward the 
observer 

Our aim in this section is to correct the positions of the par- 
ticles in the simulation according to the perturbations in the 
light cone and to obtain their "observable" coordinates. A given 
particle will be observed when it intersects the light cone of 
the observer at a certain fj. The unperturbed photon path in 
parametrized by r{ri) = i]o ^ V and constant angular coordi- 
nates that coincide with Oaiv), faiv)- The intersection will oc- 
cur when ITol[T2li 

ra + 2 ( (j>Ndri, (44) 
J no 
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daiv) = 6 + 60 = 0-2 



'^'^rf^lltZ^ (45) 
r{ij)x de 



(fiaiv) = (p + Sip = ip 



l-r(ft) 

2 / dx- 

Jo 



r{fj) - X 



'ri'ij)xsm'^ 0aifl) d(p ' 
(46) 

where 770 is the conformal time at the origin and the integrals are 
taken along the unperturbed light cone (Born approximation). 

The observed redshift of the particle is also different from the 
one that we would measure in a homogeneous universe. The 
transformation is given by conservation of energy, Eq. (O. This 
allows us to write the set of observable coordinates of the parti- 
cles as 
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0a{r])- 50, (48) 
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where V indicates the peculiar velocity projected on n. The 
terms 0Ar(O) and V{0) produced by the gravitational potential 
and velocity of the observer contribute to the monopole and 
dipole anisotropics making these terms in practice not useful as 
cosmological probes. 



VII. SUMMARY 



We have given a dictionary for how to use the outputs of nu- 
merical simulations run using Newtonian dynamics to compute 
the clustering properties of matter even on scales comparable 
to the horizon. We have shown that as long as there is a large 
separation between the length scale at which the comoving cur- 
vature ( starts to evolve with time and the scale of the horizon, 
the output of calculations based on Newtonian dynamics can be 
used even on very large-scales provided one reinterprets the co- 
ordinates of the particles. In standard large-scale structure sim- 
ulations this separation of scales results from the fact that the 
nonlinear scale is well inside the horizon, but in general it will 
occur if all species that cluster are non-relativistic and the den- 
sity perturbations are small. We gave formulas to compute the 
coordinates of particles in the conformal Newtonian gauge given 
the outputs of a simulation and to correct their positions to ob- 
servable coordinates from the same outputs. 
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Appendix: Conservation of tlie comoving curvature 



Consider a cosmological fluid with energy-momentum tensor 



^0 



{p + p)u-^up+p5-^. 



(A.l) 



a given equation of state, p[p), and speed of sound, = 
The evolution equation for the potential is given by replacing in 
Eq. (doll. 



9N 



AttGo^Sp, (A.2) 



where SpSj = —STj ai-e the pressure fluctuations. For adiabatic 
perturbations, 5p = c^Sp, then 



6^ + 3H(l + c2)0'^ 
+ [2H' + (1 + 3cl)'H^](j)N 



= 0. 



(A.3) 



In the previous equation, Cg/His the size of the sound horizon. 
Consequently, the term of order c^fc^, when compared to terms 
of order ^ TH? or ^ c^TL^, is only relevant when the typical size 
of the perturbation is smaller than the sound horizon. 

Long wavelength solutions to Eq. (IA.3b . characterized by 
kcsf] ^ 1, are easier to address in terms of a conserved quantity, 
the comoving curvature (, given by 1I22I1 



. 2 7r^ 

^ 3 1+uj 



ON- 



(A.4) 



Following I2M, to prove that the comoving curvature is con- 
served we define a new variable 



u = exp 



^ 1(1 + cDndr^ 



(A.5) 



The evolution equation obtained for u from Eq. ( IA.3I I is then 

6" 

u" - clV^u - — u = 0, (A.6) 



where 6 = ^[§(1 — ^)] In the case of long wavelength 
perturbations, the solution to the evolution equation is 



u = Aie + A2e I § 



(A.7) 



where Ai and A2 are constants. It can be shown that 



G)' 



(A.8) 



reduces to the same expression as Eq. (IA.4I ) and remains con- 
stant outside of the sound horizon. 
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Connection between Newtonian simulations and general relativity 
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On large-scales, comparable to the horizon, the observable clustering properties of galaxies are affected by 
various general relativistic effects. To calculate these effects one needs to consistently solve for the metric, 
densities and velocities in a specific coordinate system or gauge. The method of choice for simulating large- 
scale structure is numerical N-body simulations which are performed in the Newtonian limit. Even though one 
might worry that the use of the Newtonian approximation would make it impossible to use these simulations to 
compute properties on very large-scales, we show that the simulations are still solving the dynamics correctly 
even for long modes and we give formulas to obtain the position of particles in the conformal Newtonian gauge 
given the positions computed in the simulation. We also give formulas to convert from the output coordinates 
of N-body simulations to the observable coordinates of the particles. 

PACS numbers: 98.80.-k,98.80.Jk,98.62.Py,98.65.-r 



I. INTRODUCTION 



o 

(N 

CD 

o 

■ The study of the fluctuations in the distribution of matter in the 
^> Universe and its evolution through cosmic history has become 
r \one of the major tools in cosmology. The properties and time 
dvolution of the large-scale structure depend on the cosmologi- 
cal parameters and on the initial conditions for the hot big bang. 
OMany of the parameters of the currently favored cosmological 
Q model have been determined by matching the observed proper- 
?— I ties of the distribution of mass through cosmic history with the 
c/) model calculations. 

^ ^ ^ Galaxies serve as tracers of the underlying matter distribution. 

Significant efforts have been made to understand their connec- 
(T^tion iQl-Q] and to generate estimates for cosmological parame- 
J>ters from the recovered matter power spectrum |5]. In the last 
decades, redshift surveys such as the Sloan Digital Sky Survey 
•^(SDSS) (@] and flie Two-degree Field Galaxy Redshift Survey 
XJ(2dFGRS) fh have resulted in detailed maps of the large-scale 
• distribution of galaxies across very large volumes. The future 
promises even larger surveys as a result of efforts to improve 
measures of the so-called baryon acoustic oscillation signal in 
I the clustering of matter to further constrain the properties of the 
• .dark energy |l8l[9[]. Surveys are beginning to probe scales com- 
parable to the horizon at the redshift of the galaxies being ob- 
k> served. 

^ Until recently, predictions for observables in galaxy surveys 
5^ had been done entirely in the Newtonian limit. Matsubara iflOll 
included gravitational lensing effects on the correlation func- 
tions of galaxies and quasars as applied to SDSS. More recently, 
work by Yoo et al. 11 1.1 . Yoo LlZl made a more detailed treat- 



ment of general relativistic effects. These become relevant as 
the scales probed by the survey approach the horizon scale. The 
overdensity in the galaxy distribution Jobs is given by 



^obs 



h {5m - 



Z5z 
d_ 

dz 

1 + zdH 



A + 2D + {v' 



Sz-2 



H dz 



5z + 2 



l + z 
Hr 
Sr 



6z 



- 5p 6Vl 



E, 



e e' 



(1) 



where H is the Hubble constant; 5z, Sr and SVl are the fluc- 
tuations in the redshift, distance along the line of sight, and lu- 
minosity distance relative to the unperturbed universe; k is the 
lensing convergence; p gives the slope of the galaxy luminos- 
ity function; b is the bias; the direction of propagation of the 
photon; and A, Bi, D and Eij are metric components: 



ds' = 



a' [(14 



2 A) dif -2a^Bi dr] dx' 



(2) 



2D)g,^ 



2Eii\ dx'dx' 



with gij the metric tensor for three-space in a homogeneous 
universe and ?/ = J dt/a{t) is the conformal time in terms of 
the scale factor. These formulas exhibit many of the relativis- 
tic effects that are common in calculation of the anisotropics in 
the cosmic microwave background (CMB). For example the ob- 
served redshift of a source is given by 



1 + Zobs 



- 1 



B,) e'-A 



(3) 



dr[{A' -D')-{B,y+El^)e 
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where the prime indicates the derivative with respect to confor- 
mal time, the vertical bar is the covariant derivative with respect 
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to gij, Ts — r{zs) is the comoving line of sight distance to the 
source galaxies at Zg, u^e* is the line of sight peculiar veloc- 
ity, and flo and are the values of the scale factor at the time 
of observation and light-emission respectively. The first square 
bracket represents the redshift-space distortion by peculiar ve- 
locities, frame dragging, and gravitational redshift, respectively. 
The first round bracket in the integral also represents the gravi- 
tational redshift, arising from the net difference in gravitational 
potential due to its time evolution for the duration of photon 
propagation, and this effect is referred to as the integrated Sachs- 
Wolfe (ISW) effect in the CMB literature. The last terms in the 
integral represent the tidal effect from the frame dragging and 
the ISW effect from the time evolution of the primordial gravity 
waves. 

The complete set of formulas needed to predict the observed 
clustering properties of galaxies on very large-scales can be 
found in Yoo et al. [11]. It is clear that the calculation requires 
consistently solving the general relativity dynamics in a partic- 
ular coordinate system or gauge. On the other hand N-body 
simulations are the method of choice to compute predictions for 
the large-scale distribution of galaxies but these simulations are 
done in the Newtonian limit. It is appropriate to ask how the 
output of simulations can be used to compute the different terms 
in Eq. ([U and even whether this can be done at all given that the 
simulations are run using Newtonian dynamics. 

The drive on the observational side to map larger and larger 
volumes of the Universe and the exponential increase in com- 
puter power have also resulted in numerical simulations of 
ever increasing size. Typical cosmological simulations evolve 
the particles starting at z ~ 100, when the size of the hori- 
zon is ^ l.SGpc. Box sizes vary and can be as large as 
^ 0.5 — 3Gpc comoving and the number of particles is of order 
^ 10^ — 10^*^. Examples of some of the biggest simulations to 
date are the Millennium Simulation 1 13] run in a box of comov- 
ing size 500/i^^Mpc, the Marenostrum Numerical Cosmology 
Project 1 14] and the Hubble Volume project 11151] run in a box 
of 3000 X 3000 X 30 Mpc. Further examples are found in Col- 
berg et al. iH], Park et al. |[l3l and recently in Cai et al. ifisll . 
Some of these simulations are started at an initial time when the 
horizon actually lies within the box. Clearly, we need a way to 
match cosmological N-body simulations with the general rela- 
tivistic variables in Eq. ([TJ. 

In addition to asking how to use the outputs of numerical sim- 
ulations to compute the various terms in Eq. ([T]i one may won- 
der if numerical simulations are solving the correct dynamical 
equations. We might suspect that the Newtonian simulations are 
working in the so-called conformal Newtonian gauge, in which 
the line element is given by 

ds^ = a^(?7)[-(l + '2.(t)N)drf + (1 - 2(j)N)5^jdx'dx^] (4) 

in the absence of anisotropic stress and where 4)^ coincides with 



the Newtonian potential only on small scales. In fact, the ana- 
logue of the Poisson equation in the conformal Newtonian gauge 
reads: 

V'^(j)N -in{(l)'N+n(l)N)=^T^Ga^5pN (5) 

and thus differs from the standard Poisson equation on large- 
scales. Here, T-L — a' /ais the conformal Hubble parameter. 

From Eq. (|5]l, it might appear that simulations are not solving 
correctly for the gravitational potential for scales comparable to 
or larger than the horizon. Previous work on this subject has fo- 
cused on comparing the general relativity equations to the New- 
tonian equations up to some given order in perturbation theory 
ifigll . We will show in this work that a more direct approach is 
possible. We will analyze the situation in detail and conclude 
that N-body simulations are solving for the potential correctly 
but that the location of the particles needs to be corrected if they 
are to be interpreted as the particle coordinates in the conformal 
Newtonian gauge. Finally we will give formulas to recover ob- 
servable coordinates directly in terms of the output of N-body 
simulations. 



II. EVOLUTION EQUATIONS 

As structure in the Universe develops the density contrast be- 
comes larger and larger exceeding unity at the so-called nonlin- 
ear scale. Properly modeling this process on small scales, of 
order the nonlinear scale or smaller requires numerical simula- 
tions. However, because the primordial curvature fluctuations, 
the seeds for structure formation, are so small, the nonlinear 
scale is significantly smaller than the horizon. As a result per- 
turbations in the space time remain very small, of order 10^^ or 
smaller'. 

Thus, to study structure formation we need only consider 
small perturbations to the Friedmann-Robertson- Walker metric 
and we can stay at linear order on those perturbations. In this pa- 
per we choose to work in the conformal Newtonian gauge with 
the line element given by: 

ds^ = a2(?/)[-(l + 2ijN)d7f + (1 - 2(l)N)5^jdx'dx^] (6) 

where ipM represents the Newtonian potential and (f)^, the New- 
tonian curvature. 



' On sufficiently small scales baryons can collapse to fomi relativistic objects 
such as neutron stars or black holes around which the space time metric fluc- 
tuations are large. This has negligible effects on the length scales considered 
in this paper 
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We stress that we are assuming that metric perturbations are 
small but we are not treating the density perturbations using per- 
turbation theory. The structure formation process also results in 
peculiar velocities for the particles. Because the nonlinear scale 
is well inside the horizon, these peculiar velocities are small, 
much smaller than the speed of light. In fact at the nonlinear 
scale peculiar velocities are of order the Hubble velocity for 
points separated by a distance of order the nonlinear scale. As 
a result the kinetic energies of particles do not source gravity in 
an appreciable way. 

Note that we are considering perturbations around the FRW 
metric so at lowest order the source for the gravitational poten- 
tials in 5p as opposed to the full p. The kinetic energy correc- 
tions are of order pv^. It is still true that pv^ <^ 6p on every 
scale. Thus it is safe to ignore the peculiar motions as a source 
of gravity. These terms are of course also neglected in numerical 
simulations run using the Newtonian approximation, but this is 
a negligible source of error. Including these terms is necessary 
if one wants to study the back-reaction of cosmological pertur- 
bations on to the expansion of the Universe [20], but we are not 
interested in this problem here. 

In the standard Newtonian approximation terms of order p(j) 
are also dropped as sources of gravity. This requires a little bit 
more thought in our case. Again at lowest order the source of 
gravity is 6p but this is no longer much larger than p(p on suffi- 
ciently large-scales. Thus we need to keep this term. However 
it is only p(f) that needs to be kept as of course Sp 3> Spcf) on all 
scales. We will now summarize the evolution equations under 
these approximations. 



A. Einstein equations 

In the conformal Newtonian gauge the Einstein equations 

Gf_i,y — SnGT^f, are reduced to: 



where ^ indicates derivatives with respect to the i coordinate 
and the background Friedmann equation for a flat Universe with 
cosmological constant A has been subtracted. 



— = -A.GT,+-, 



(11) 



B. Application to nonrelativistic particles 

As we mentioned when studying structure formation we are 
primarily interested in nonrelativistic matter Equation ( fTOb for 
i ^ i implies that the anisotropic stress is of order pv'^ and thus 
negligible in our approximation. As a result (pN = ipN- The 
gravitational potential satisfies: 



3n( 



HM = -ATiGa\T^ - f-O). (12) 



The energy-momentum tensor for a set of cold dark matter par- 
ticles with mass ma [21] is given by: 



T^" = i-g) 



,-1/2 "a"a 

g] ' y m, 



5d{x -Xa) 



(13) 



where is the comoving velocity of the particles, dx'^'' /drj and 
g is the determinant of the metric. At linear order in the metric 
and to second order in the three-velocities, u*, this is given by 



^0 = a-i (1 - + vl/2) and < = 



The 00 compo- 



nent of the stress-energy tensor is related to the density of parti- 
cles and because metric perturbations are small, we can expand 
in powers of and remain at linear order. To order v^, the 00 
component is 

Tl^^-a-^'Y^mail + S^N+vDdDix^-Xam- (14) 



V^cbN - m{4>'N + n^N) - -47rGa2(T° - f^), (7) 



b'^+n^Nh^^TiGa^Tl 



(8) 



As discussed we will neglect the term but need to keep the 
4)N as (fiMP is not negligible on large-scales. We obtain 

T° = -a-^(l + 30jv) "^aSoix: - Xa{v)) (15) 



AtiGo^ 



1 



hj^'^{4>N - i'N)] 



(9) 



(10) 



where only the 0jvP piece of the term proportional to (j)N ever 
makes any difference. Replacing in Eq. (fT2l i. 



47rGaV(l + 3(^Ar) + — 



(16) 



where p{x, rf) — a ^ niaSDix — Xa{ri)) is the density ob- 
tained by naively counting particles in cells at each time step. 
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Given the positions of the particles, Eq. ( fTSI l can be solved to 
obtain the Newtonian potential. 

The positions and velocities of the particles are advanced us- 
ing the geodesic equation. 



drj 



(17) 



Notice that the term 3(j)'j^dxa/dri is always negligible. 

m. INITIAL CONDITIONS 

In addition to the evolution equations we need to find the ini- 
tial conditions. This can be done at early enough time using 
linear theory. We define a growth function for the potential such 
that in the linear regime 



(t>N = b^{vi)(j)' 



N- 



(18) 



Given that the spatial components of the stress tensor for the 
dark matter are negligibly small, the gravitational potential in 
the matter era satisfies 



N 



0. 



(19) 



The solutions to this equation are well-known f^^], giving a 
constant and a decaying solution for the potential in the linear 
regime, which in terms of Eq. ( fTSl ) is 



9i 



(20) 



where Ci and C2 are constants. The decaying mode is abso- 
lutely negligible at the times of interest and without loss of gen- 
erality we choose Ci — 1. 

We can replace the right-hand side of the geodesic equation 
by the solution of the potential in the linear regime. 



^dXa ^ _ 

dr] 



in 
N- 



(21) 



In a matter-dominated regime the homogeneous solutions to Eq. 
(ETJ are given by a constant vector and a decaying solution. 



Xh Bir]-^ + B2. 



(22) 



For 



the particular solution, we choose an ansatz 
— bs{r]) ijMxin), as it is usually done in the Zel'dovich 
approximation |23]. The labeling of 65(77) as such will become 
clear by the end of this section. The equation to be solved is. 



N- 



(23) 



The right-hand side of the last equation is independent of time, 
which implies that h'^ + T-ib'g — constant. The actual value 
of the constant is arbitrary, since it can be absorbed in ipi{xin)- 
Thus, bs oc f]^. (Adding a constant to this solution would not 
modify the subsequent steps of our paper and is only linked to 
the choice of initial time.) If we equate the factors that depend 
on the coordinates, ^pl{xin) oc —V(j)^j!}. 

To give the complete solution for the position of the particles 
we separate the constant vector B2 in two components, B2 = 
Xin + 6xin{xin), whcrc Xin are the positions of the particles if 
they started out distributed uniformly in a mesh. We discard the 
decaying term and give the position of the particles as a function 
of time in the linear regime. 



Xa — ^•^in 

(fj„) + bs{rj)%IJi{xin). 



(24) 



To find the value of Sxin{xin), we resort to the Poisson equation 
(fTSI l. We can evolve p{x, rf) by means of a transformation of 
coordinates from the initial particle density 



p(x, ri) 



P 



(25) 



where 



dx 



is the Jacobian of the transformation and p is the 



initial uniform background density. The transformation of co- 
ordinates is given by Eq. ( |24] |. where 6x, bs{r]) and -01 are un- 
knowns. Since the perturbations are initially small, the density 
evolves as 



P 



(l-V-(52;„-&5(?7)V-Vi). 



(26) 



We can define contributions to the density perturbation, 5, as 
related to the displacement fields by 



-V • 5xi 



'bs(rj)V ■ -01 



(27) 



(28) 



The physical meaning of bs{ri) now becomes clear, as it is iden- 
tified with the growth function of density perturbations. Indeed, 
as we expect in the matter-dominated regime, we obtained that 
65(77) oc 0(77) oc ?7^. 

To determine we replace Eqns. (l26b . (l27b and (|28]) in Eq. 
(fTSI l to obtain: 



S^). (29) 



At the initial time the first term in the left-hand side cancels with 
the rightmost term in the right-hand side of the previous equa- 
tion. We can solve for (5"' from the remaining terms. 



N- 



(30) 
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In terms of their Fourier components, while (5™ oc (/i™ j,, the 
Zel'dovich component dependency is oc fc^H^^i/)^^ Well 
inside the horizon, in the limit kr] ^ 1, (5™ becomes negligible 
as compared to , the Newtonian density perturbation. When 
kr] 1 then ^ and as anticipated we cannot neglect 
this term. 

Finally, we can obtain 6xin by inverting Eq. (IZTl i in Fourier 
space. 



(Pk ik 



(31) 



In summary, to perform a cosmological simulation we need to 
evolve the position and velocities of particles using Eq. (fTTI ) and 
compute the potential using Eqs. Q and ( fTSl l. At the initial time 
the cold dark matter particles need to be displaced by an amount 
given in Eq. dSTT i. 



IV. COMPARISON TO NEWTONIAN COSMOLOGICAL 
SIMULATIONS 

Now that we have a consistent set of equations to solve we 
can compare them to those used in cosmological simulations to 
determine whether these simulations can be used to study very 
long wavelength modes or if they require some change. 

Cosmological Newtonian simulations solve for the potential 
by means of the Poisson equation 



(32) 



and move particles according to Newton's law expressed in co- 
moving coordinates. 



(Fx a 



di] 



(33) 



This evolution equation is identical to the geodesic equation ( fTTI ) 
if the gravitational potential were computed correctly (as men- 
tioned before the term oc S^at' is negligible for modes both large 
and small compared to Hubble). Thus if the gravitational poten- 
tial is correct the particle positions are updated properly. 

It is important to determine if there are corrections to the grav- 
itational potential that become important on large-scales. The 
density that sources the Poisson equation in simulations is di- 
rectly computed by counting particles in cells, 

Psim ix,r/) ^ Msim ^ a'^^maSoix - Xail])). (34) 

a 

Even if the particle positions had been computed correctly, this 
"simulation density" differs from the density in the conformal 
Newtonian gauge by a factor (1 + S^at). 



Finally in standard cosmological simulation the particles are 
initially displaced making use of the Zel'dovich approximation, 
which in the Newtonian case takes the form. 



Xa = Xin + bs{r])'ljji{xin). 



(35) 



This differs from the displacements we calculated in the previous 
section - it is missing the Sxin- 

Thus at first sight it appears that the gravitational potential is 
not computed using the correct equation, that the density con- 
trast is missing a term and that the initial displacement of the 
particles is incorrect. We will now show that in fact all these dif- 
ferent "missing terms" cancel each other so that the gravitational 
potential is computed correctly. As a result, particle positions are 
also updated correctly. 

Let us look at the situation more carefully. For completeness 
let us also include a cosmological constant and start by restrict- 
ing ourselves to the linear regime as in any event the effects we 
are considering are only relevant on very large-scales. The rela- 
tivistic Poisson equation reads 



(36) 



where we have introduced the equation of state parameter co — 
p/p and we have written —V • Sxin = S"^ = SCin- Notice that 
3/2H^(l + u) = ATrGpdm with pdm the mean density for the 
matter. The terms in brackets on the right-hand side correspond 
to the density contrast calculated in the Newtonian simulations 
{5^) and the two missing corrections, the one proportional to 
and the one coming from the missing initial displacements (Cin)- 
It is useful to consider the comoving curvature defined as 



c 



27r 
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1 



(37) 



It is well known that this comoving curvature remains constant in 
time on large-scales. For completeness we spell out the deriva- 
tion in the Appendix. Note that C,in defined above is nothing 
other than the initial value for this variable In the case of 
a perfect fluid C is constant on all scales larger than the sound 
horizon, or fc^c^ << TH? . In general what plays the role of 
is just the typical velocity dispersion that relates the magnitude 
of the spatial components of the energy-momentum tensor to the 
density. In our case, it is the velocity dispersion of the dark mat- 
ter particles induced by the growth of perturbations and thus 
remains constant all the way to the non linear scale. 

The relativistic corrections to the Newtonian Poisson equa- 
tion, V^0Ar — 3/2H^(l + oj)5^ , are evidenced from subtract- 
ing this expression from Eq. (l36T l. Notice that the difference 
between the relativistic terms in the right- and left-hand sides of 
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the Poisson equation is nothing other than: 

3H(0^ + n4>N) + ^n^i + u){4>N - Cm) = (38) 

Thus the additional terms in the relativistic Poisson equation 
cancel each other for modes larger than the nonlinear scale. Fur- 
thermore, because the nonlinear scale is well inside the horizon 
once the cancellation begins to fail the additional terms are very 
small, of order the ratio of the nonlinear scale to the horizon 
squared. Thus the gravitational potential in Newtonian simula- 
tions coincides with the one in the conformal Newtonian gauge 
even on very large-scales and thus the dark matter particles are 
being displaced correctly. 

We have showed that Newtonian N-body simulations calcu- 
late the correct gravitational potential and displace the particles 
correctly. The coordinates of the particles however are not the 
coordinates in the conformal Newtonian gauge as they are miss- 
ing the initial displacement. Thus the dictionary between con- 
formal Newtonian gauge variables and numerical simulations is: 

</>Af = 0sim, (39) 
m = Vsim, (40) 
XN = Xsim + Sxin- (41) 

Notice that the reason why the scheme worked was that the 
"sound horizon" scale, the scale out to which C was constant in 
time, was well inside the horizon. This scale is nothing other 
than the scale that dark matter particles can move since the time 
of the Big Bang as a result of the peculiar velocities they have. 
Because the dark matter particles are non-relativistic this dis- 
tance is well inside the horizon and the mistakes are negligible. 
Of course a simulation based on Newtonian physics could not 
work if particles are moving at an appreciable fraction of the 
speed of light. In general if there is an additional relativistic 
component the Newtonian simulations would not be computing 
things properly. The cosmological constant did not cause any 
problem because even though it is in some sense relativistic it 
is homogeneous so it does not contribute to the perturbation of 
the stress tensor. Thus as long as one is modeling nonrelativis- 
tic components or a relativistic component that does not cluster 
one is safe using the Newtonian approximation. Such a New- 
tonian approximation would not work, for example, during the 
radiation era where the effective sound speed of the dominant 
component of the energy density is very close to the speed of 
light. 



V. THE COMOVING GAUGE 

For completeness we will now show that Eq. ( fT2] l can be writ- 
ten making apparent gauge transformations between the comov- 
ing and the conformal Newtonian gauge. 

In Eq. ( l36b . the factor -^S^ is the density perturbation in 
Newtonian simulations. The expression in brackets is the con- 
formal Newtonian gauge density perturbation. Indeed rearrang- 
ing the terms we have 

\7^^N-5ni(b'N + nM = (42) 

As long as ( is constant, we can recognize in the right-hand side 
the gauge transformation for the density perturbation between 
the comoving and the conformal Newtonian gauge, 

SpN^5pc + S{^N-C)pi^+^)- (43) 

This is a well-known relation ll24ll and in this context it implies 
that the density perturbation that Newtonian simulations are ob- 
taining is the one in the comoving gauge. 

VI. OBSERVABLE COORDINATES 

In an inhomogeneous universe, the observed positions of the 
particles in the simulation are modified due to effects such as 
the Sachs-Wolfe effect, gravitational lensing, and peculiar ve- 
locities. The net result is that photons from a source follow a 
path that is perturbed with respect to the light cone of an ob- 
server in a homogeneous universe. Consider a comoving ob- 
server in an inhomogeneous universe with a velocity given by 
— {{1 — iI>m) / cLtV^ / a). The direction of observation is n, 
defined by {6, (p) in spherical coordinates, but due to the per- 
turbations to the photon path, the direction toward the point 
reached by the photon geodesic is actually s, corresponding to 
{9 + S9,(p + 6(p). We now take a particle with coordinates Xa{ri) 
(already taking into account the correction Sxin) and we want to 
know where it crosses the path of the photons going toward the 
observer 

Our aim in this section is to correct the positions of the par- 
ticles in the simulation according to the perturbations in the 
light cone and to obtain their "observable" coordinates. A given 
particle will be observed when it intersects the light cone of 
the observer at a certain fj. The unperturbed photon path in 
parametrized by r{ri) = i]o ^ V and constant angular coordi- 
nates that coincide with Oaiv), faiv)- The intersection will oc- 
cur when ITol[T2li 

ra + 2 ( (j>Ndri, (44) 
J no 
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daiv) = 6 + 60 = 0-2 



'^'^rf^lltZ^ (45) 
r{ij)x de 



(fiaiv) = (p + Sip = ip 



l-r(ft) 

2 / dx- 

Jo 



r{fj) - X 



'ri'ij)xsm'^ 0aifl) d(p ' 
(46) 

where 770 is the conformal time at the origin and the integrals are 
taken along the unperturbed light cone (Born approximation). 

The observed redshift of the particle is also different from the 
one that we would measure in a homogeneous universe. The 
transformation is given by conservation of energy, Eq. (O. This 
allows us to write the set of observable coordinates of the parti- 
cles as 
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Zo,. = '-^[l + Vis) 
a{ri) 



Oobs 
P>obs 



V{Q)~<t>N{s)+4>Nm (47) 

r(jj) 



1, 



0a{r])- 50, (48) 
iPaiv)-6p. (49) 



where V indicates the peculiar velocity projected on n. The 
terms 0Ar(O) and V{0) produced by the gravitational potential 
and velocity of the observer contribute to the monopole and 
dipole anisotropics making these terms in practice not useful as 
cosmological probes. 



VII. SUMMARY 



We have given a dictionary for how to use the outputs of nu- 
merical simulations run using Newtonian dynamics to compute 
the clustering properties of matter even on scales comparable 
to the horizon. We have shown that as long as there is a large 
separation between the length scale at which the comoving cur- 
vature ( starts to evolve with time and the scale of the horizon, 
the output of calculations based on Newtonian dynamics can be 
used even on very large-scales provided one reinterprets the co- 
ordinates of the particles. In standard large-scale structure sim- 
ulations this separation of scales results from the fact that the 
nonlinear scale is well inside the horizon, but in general it will 
occur if all species that cluster are non-relativistic and the den- 
sity perturbations are small. We gave formulas to compute the 
coordinates of particles in the conformal Newtonian gauge given 
the outputs of a simulation and to correct their positions to ob- 
servable coordinates from the same outputs. 
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Appendix: Conservation of tlie comoving curvature 



Consider a cosmological fluid with energy-momentum tensor 



^0 



{p + p)u-^up+p5-^. 



(A.l) 



a given equation of state, p[p), and speed of sound, = 
The evolution equation for the potential is given by replacing in 
Eq. (doll. 
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AttGo^Sp, (A.2) 



where SpSj = —STj ai-e the pressure fluctuations. For adiabatic 
perturbations, 5p = c^Sp, then 



6^ + 3H(l + c2)0'^ 
+ [2H' + (1 + 3cl)'H^](j)N 



= 0. 



(A.3) 



In the previous equation, Cg/His the size of the sound horizon. 
Consequently, the term of order c^fc^, when compared to terms 
of order ^ TH? or ^ c^TL^, is only relevant when the typical size 
of the perturbation is smaller than the sound horizon. 

Long wavelength solutions to Eq. (IA.3b . characterized by 
kcsf] ^ 1, are easier to address in terms of a conserved quantity, 
the comoving curvature (, given by 1I22I1 



. 2 7r^ 

^ 3 1+uj 



ON- 



(A.4) 



Following I2M, to prove that the comoving curvature is con- 
served we define a new variable 



u = exp 



^ 1(1 + cDndr^ 



(A.5) 



The evolution equation obtained for u from Eq. ( IA.3I I is then 

6" 

u" - clV^u - — u = 0, (A.6) 



where 6 = ^[§(1 — ^)] In the case of long wavelength 
perturbations, the solution to the evolution equation is 



u = Aie + A2e I § 



(A.7) 



where Ai and A2 are constants. It can be shown that 



(A.8) 



reduces to the same expression as Eq. (IA.4I ) and remains con- 
stant outside of the sound horizon. 



[1] N. Kaiser, Astrophys. J. Lett. 284, L9 (1984). 

[2] M. J. Rees, Mon. Not. R. Astron. Soc. 213, 75P (1985). 

[3] 1. M. Bardeen, 1. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. 

J. 304, 15 (1986). 
[4] H. 1. Mo and S. D. M. White, Mon. Not. R. Astron. Soc. 282, 347 

(1996), arXiv:astro-ph/95 12127. 
[5] B. A. Reid, W. 1. Percival, D. 1. Eisenstein, L. Verde, D. N. 

Spergel, R. A. Skibba, N. A. Bahcall, T. Budavari, 1. A. Frieman, 

M. Fukugita, et al., Mon. Not. R. Astron. Soc. 404, 60 (2010), 

0907.1659. 

[6] D. G. York, 1. Adelman, J. E. Anderson, Ir., S. F. Anderson, 
J. Annis, N. A. Bahcall, 1. A. Bakken, R. Barkhouser, S. Bas- 
tian, E. Berman, et al., Astron. 1. 120, 1579 (2000), arXiv:astro- 
ph/0006396. 

[7] M. Colless, G. Dalton, S. Maddox, W. Sutherland, R Norberg, 
S. Cole, 1. Bland-Hawthorn, T. Bridges, R. Cannon, C. Collins, 
et al., Mon. Not. R. Astron. Soc. 328, 1039 (2001), arXiv:astro- 
ph/0106498. 

[8] D. 1. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, M. R. 
Blanton, R. C. Nichol, R. Scranton, H. Seo, M. Tegmark, 
Z. Zheng, et al., Astrophys. J. 633, 560 (2005), arXiv:astro- 
ph/0501171. 



[9] W. J. Percival, B. A. Reid, D. J. Eisenstein, N. A. Bahcall, T. Bu- 
davari, J. A. Frieman, M. Fukugita, 1. E. Gunn, Z. Ivezic, G. R. 
Knapp, et al., Mon. Not. R. Astron. Soc. 401, 2148 (2010), 
0907.1660. 

[10] T. Matsubara, Astrophys. J. Lett. 537, L77 (2000), arXiv:astro- 
ph/0004392. 

[11] J. Yoo, A. L. Fitzpatrick, and M. Zaldarriaga, Phys. Rev. D 80, 
083514 (2009), 0907.0707. 

[12] J. Yoo, Phys. Rev. D 82, 083508 (2010), 1009.3021. 

[13] V. Springel, S. D. M. White, A. lenkins, C. S. Frenk, N. Yoshida, 
L. Gao, J. Navarro, R. Thacker, D. Croton, J. Helly, et al.. Nature 
435, 629 (2005), arXiv:astro-ph/0504097. 

[14] S. Gottlober, G. Yepes, A. Khalatyan, R. Sevilla, and V. Turchani- 
nov, in The Dark Side of the Universe, edited by C. Manoz & 
G. Yepes (2006), vol. 878 of American Institute of Pliysics Con- 
ference Series, pp. 3-9, arXiv:astro-ph/0610622. 

[15] A. E. Evrard, T. J. MacFarland, H. M. P Couchman, J. M. Colberg, 
N. Yoshida, S. D. M. White, A. Jenkins, C. S. Frenk, F. R. Pearce, 
J. A. Peacock, et al., Astrophys. J. 573, 7 (2002), arXiv:astro- 
ph/0110246. 

[16] J. M. Colberg, S. D. M. White, N. Yoshida, T. J. MacFarland, 
A. Jenkins, C. S. Frenk, F. R. Pearce, A. E. Evrard, H. M. P. 



9 



Couchman, G. Efstathiou, et al., Mon. Not. R. Astron. Soc. 319, 

209 (2000), arXiv:astro-ph/0005259. 
[17] C. Park, J. Kim, and J. R. Gott, III, Astrophys. J. 633, 1 (2005), 

arXiv:astro-ph/0503584. 
[18] Y. Cai, S. Cole, A. Jenkins, and C. S. Frank, Mon. Not. R. Astron. 

Soc. 407, 201 (2010), 1003.0974. 
[19] J. Hwang and H. Noh, Mon. Not. R. Astron. Soc. 367, 1515 

(2006), arXiv:astro-pli/0507159. 
[20] D. Baumann, A. Nicolis, L. Senatore, and M. Zaldarriaga, ArXiv 



e-prints (2010), 1004.2488. 

[21] S. Weinberg, Gravitation and Cosmology: Principles and Appli- 
cations of the General Theory of Relativity (1972). 

[22] V. Mukhanov, Physical Foundations of Cosmology (2005). 

[23] Y. B. Zel'dovich, Astron. Astrophys. 5, 84 (1970). 

[24] W. Hu, ArXiv Astrophysics e-prints (2004), arXiv:astro- 
ph/0402060. 



